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Reliable noise prediction capabilities are essential to enable novel fuel efficient open ro- 
tor designs that can meet the community and cabin noise standards. Toward this end, im- 
mersed boundary methods have reached a level of maturity so that they are being frequently 
employed for specific real world applications within NASA. This paper demonstrates that 
our higher-order immersed boundary method provides the ability for aeroacoustic analysis 
of wake-dominated flow fields generated by highly complex geometries. This is the first of 
a kind aeroacoustic simulation of an open rotor propulsion system employing an immersed 
boundary method. In addition to discussing the peculiarities of applying the immersed 
boundary method to this moving boundary problem, we will provide a detailed aeroacous- 
tic analysis of the noise generation mechanisms encountered in the open rotor flow. The 
simulation data is compared to available experimental data and other computational re- 
sults employing more conventional CFD methods. The noise generation mechanisms are 
analyzed employing spectral analysis, proper orthogonal decomposition and the causality 
method. 


I. Introduction 


There has been a renewed interest in pursuing contra-rotating open rotor (CROR) propulsion technology 
due to the large potential of significantly reducing fuel consumption. A main concern for the design of such 
systems is that they must meet community noise and cabin noise standards. Hence, providing reliable noise 
prediction capabilities for contra-rotating open rotors is essential in the design of low-noise rotor propulsion 
systems. In recent years, NASA initiated several research efforts that have assessed the capabilities on noise 
source prediction and analysis of contra-rotating open rotors. Several low to high fidelity simulation tools 
were used to provide noise predictions and their results were compared to experimental data.b° One of the 
key challenges is to conduct efficient high-fidelity simulations of this relatively complex geometry including 
moving boundaries. 

While providing a detailed acoustic analysis of the noise generation mechanisms, a key contribution of 
this paper is to conduct this analysis by employing a higher-order accurate immersed boundary method. ‘This 
is the first ever application of an immersed boundary method for noise prediction of an open rotor system. 
It will be demonstrated that this method provides a reliable noise prediction capability. Immersed Boundary 
Methods (IBMs) have been developed for many years and have appeared in various forms since they were first 
introduced by Peskin®’’ (see for example Goldstein et al.,° LeVeque and Li,? Wiegmann and Bube,’? Linnick 
and Fasel,44 Johansen and Collela,!? Mittal and Iaccarino,'? Zhong,'* Duan et al.!° and many others). These 
methods were first introduced as a nontraditional approach for numerically solving initial /boundary-value 
problems for complex geometries on Cartesian meshes and have matured to become increasingly important 
for a wide range of applications.‘° !® One of the key advantages of immersed boundary methods is that the 
computational meshes can be automatically be generated starting from a water tight surface triangulation 
independent of the complexity of the geometry. The grid-generation process for body-fitted structured or 
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unstructured grids for complex geometries is generally very laborious because the process aims at generating 
well behaved grids that have sufficient local resolution while minimizing the total number of required grid 
points. Except for simple geometries, these conflicting requirements can lead to deterioration in grid quality, 
thereby negatively impacting the accuracy and convergence properties of the flow solver. For highly complex 
geometries, which are relevant for many fields of science and engineering, the process of generating a high- 
quality grid is extremely time consuming. For flows involving moving and deforming boundaries, IBMs also 
provide clear advantages over classical body-fitted grid approaches. Simulating such flows on body-fitted 
erids requires generating a new grid at each time step and a procedure to project the solution onto this new 
erid (Tezduyar!’). These two characteristics (grid deformation and projection of solution) for body-fitted 
erids associated with simulating flows with moving and deforming boundaries may negatively impact the 
accuracy, robustness, and computational cost of the numerical solution method. Particularly, in cases where 
the boundary exhibits large motions, body-conformal grid strategies cause immense difficulties in the solution 
procedure. This is where an immersed boundary approach with its fully Eulerian approach provides great 
advantages. IBMs provide a much more convenient way of including the body motion and deformation by 
using a stationary non-deforming Cartesian grid. While these methods simplify the grid generation process, 
a detailed mathematical understanding of the IBM is necessary to avoid a negative impact of the boundary 
treatment on the robustness, the convergence behavior and the accuracy of the numerical scheme. 

Immersed boundary methods also posses well known shortcomings that have limited their applicability 
to wider range of applications. Most IBMs are only lower-order accurate and suffer from robustness issues 
when extended to higher-order. The near wall accuracy is, however, critical when simulating laminar to 
turbulent transition scenarios for wall-bounded flows to capture relevant instability mechanisms. Higher- 
order immersed boundary methods have been presented in a number of research studies by Linnick and 
Fasel,'! Zhong,'* Duan et al.,!° Brehm and Fasel,?? Brehm et al.,24 and others. The current method 
employed within NASA’s LAVA CFD solver!® is based on the higher-order immersed boundary method 
discussed in Brehm et al.?! 

The current abstract provides a rough idea of what can be expected in the final paper. The abstract 
proceeds as follows: section II introduces the numerical scheme and immersed boundary method used to 
solve the compressible Navier-Stokes equations. Section III presents the geometry used in the simulation 
and provides an overview of the simulation setup. Section IV provides an initial unsteady flow visualization 
displaying some of the complicated flow features. Section V provides a preliminary comparison of the 
present CFD data obtained on a coarse grid with available experimental data acquired in NASA’s 8-foot x 
6-foot wind tunnel. A preliminary acoustic analysis is provided in section VI employing proper orthogonal 
decomposition and spectral analysis of the pressure field in the acoustic near field. A brief description of 
the causality method is also provided in section VI that will be used to link the unsteady flow features in 
the noise source regions to the acoustic signal in the far-field. Finally, some preliminary conclusions and an 
outlook to the final paper is provided in section VII. 


II. Numerical Methods 


II.A. Interior and Temporal Discretizations 


The compressible Navier-Stokes equations considering an ideal, Newtonian, non-reactive gas are 


Op 
——— . = 1 
ay + V (pu) = 9, (1a) 
O 
“5 + V: (puu+ pd — x) =0, (1b) 
OE 
5 + V-(Butu-(pi-z)-KVT)=0, with oo 
ee aes ee 
ea ee u and p= opRT. (1d) 


In the above equations p is the density, u is the velocity vector, 6 is the unit tensor, p is the static pressure, 
T is the temperature, F is the total energy, y is the specific heat ratio, R is the gas constant, and « is the 
thermal conductivity. ‘The viscous stress tensor is 


r= 1(28) +(9-3u) (Vw) 6, (2) 
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The dynamic viscosity is given by p, @ is the bulk viscosity, and the strain rate tensor S' is defined as 
s=1 (Vu i (Vu)"). 


Equations (1)-(2) can be re-written in vector form as follows 


a +v-(F-#,) =0, with Ff = (F',F’, F*) and FSF) FF) (3) 
= 


The conservative variable vector W = (p, pu, pv, pw, pH — p)* and the inviscid fluxes are 


ieee | a | | 


pu aeap puuU Ppwu 
F' = PUuv _F? = pu? +p , and F? = PWV : 
puw puvw pw" rp 
pull pvt pw 


where the total enthalpy is H =h+wu-wu/2. The viscous fluxes are 


0 0 0 

Tox Tyx Tzx 

F,= Tay ee Tyy , and FY = Tzy 
Txz Tyz Tez 

ic fp fi; 


The interior discretization of the convective terms follows the approach presented in Laible and Fasel.?? The 
forward-flux, F*, and the backward-flux, F~, are discretized with n*”-order accurate centered upwind finite 
differences as proposed by Zhong.?° As in Laible and Fasel,?” a parameter a is used to determine the degree 
of upwinding for the grid-centered upwind differences 


Oo i+Ns aes aN a hp 
ae | = Ds ChePK = OAL SNi=1 ? (4) 
- k=i1-N, t 





where c; is the finite difference stencil coefficient, Ax is the averaged grid spacing over the spatial extent of 
the stencil and N, determines the number of grid points in the stencil. Note that a central difference stencil 
can be achieved by choosing a=0. In this case a higher-order filter based on work by Visbal and Gaitonde”* is 
used to regularize the solution. In accordance to the chosen discretization in Laible and Fasel,?? the n‘”-order 
accurate upwind stencils are calculated with n = 9,7,5,3,1 and a = —1500, 72, —12,0.6,—1. In the current 
paper, it is assumed that the numerical fluxes are approximated with flux vector splitting schemes, i.e., Lax 
Friedrichs, van Leer, AUSMPW4-, etc. The compressible Navier-Stokes equations are advanced in time with 
an explicit fourth-order accurate Runge-Kutta scheme based on the Shu-Osher formulation.?° For accurately 
capturing flow discontinuities higher-order shock capturing schemes”® are also available. The computational 
domain is decomposed into block-structured sub-domains supporting adaptive-mesh refinement as well as 
generalized curvilinear coordinates. In order to maintain higher-order accuracy across coarse-fine block 
interfaces, up-to fourth-order accurate prolongation and restriction operators are available. 


II.B. Immersed Boundary Method 


In this section, the theoretical basis of the IIM presented in Brehm and Fasel?°? is extended to the com- 


pressible Navier-Stokes equations. Many immersed methods have been developed in the past; however, in 
the derivation of these schemes, usually only the order of the local truncation error or accuracy of the nu- 
merical scheme has been considered. A posteriori, the numerical stability of these schemes is commonly 
demonstrated (in a global sense) by considering a number of different test-problems, and in a few cases, an 
additional global matrix stability analysis is employed (see for example Zhong!*). The basis of the IIM by 
Brehm and Fasel?? is that the stencil coefficients are locally optimized in order to improve the stability of 
the scheme. 

The notion of developing a strategy for improving the spectral properties of immersed methods originated 
from the idea that the stability of the numerical scheme can be formulated as an N-dimensional optimization 
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problem, where N represents the number of irregular grid points. Instead of solely considering the local 
truncation error in the derivation of the finite-difference grid stencil at an irregular grid point, an over- 
determined system of equations is constructed to determine the stencil coefficients. We employed an over- 
determined system of equations because it allows one to select free parameters, which can be tuned towards 
optimizing the local stability properties of the finite-difference stencil. Thus, assuming one free parameter 
per irregular grid point an N-dimensional optimization problem for N irregular grid points is obtained. 

The objective function of the optimization problem must be carefully chosen and depends on the nature 
of the particular partial-differential equation of interest. For advection-diffusion-type equations, the spectral 
radius of the update matrix may be considered an appropriate objective function. In the current paper, we 
refer to the update matrix as the matrix which updates the solution from time ft, to the next time-step, 
tnii. The update matrix, therefore, contains information about the spatial and temporal discretization 
of the numerical scheme. ‘The spectral radius is well-suited in this situation, because its actual value is 
associated with the stability of the numerical scheme and it can be used to formulate a necessary condition 
for the stability of the update matrix. For non-linear systems of equations, a linearization is necessary prior 
to applying the outlined linear stability analysis concepts. Note that the immersed boundary method used 
in this paper solely uses linear stability theory and, therefore, neglects possible non-linear effects on the 
stability. 

In principal, it is possible to solve a coupled N-dimensional optimization problem for N irregular grid 
points. From a practical point of view, however, it is not very efficient or desirable to solve such a large system 
of equations, especially for FSI problems where the grid topology is time dependent, thereby requiring this 
procedure to be repeated at each time-step. Therefore, in the current approach, we aimed at isolating the 
boundary stencils from the remainder of the computational domain so that the optimization problem can be 
formulated for each irregular grid-point separately. Instead of solving a global N-dimensional optimization 
problem we solve 1D optimization problems at each of the N irregular grid point individually. This approach 
ereatly reduces the computational effort needed to determine the “optimal” stencil coefficients. We refer to 
the aspect of turning the global N-dimensional problem into N one-dimensional problems as “localization” 
or “localization assumption”. In Brehm and Fasel,”° it was demonstrated, both numerically and analytically, 
that for advection-diffusion type equations, the localization assumption appears to be valid as long as the von 
Neumann number does not reach a limiting value. For a more detailed discussion on the original approach, 
the interested reader is referred to Brehm and Fasel.2° In the following discussion, the basic steps for 
extending the ITM for the compressible Navier-Stokes equations are outlined. 

The main difference for the application of the IIM to the compressible Navier-Stokes equations instead 
of the incompressible Navier-Stokes equations (solved by utilizing the approximate projection method) is 
the use of numerical fluxes. Flux vector splitting schemes, i.e. van Leer,?” Rusanov, and AUSMPW4+, are 
adopted here. In the flux vector splitting approach, the convective flux is divided into a forward-flux and 
a backward-flux, F = F~ + F*. The flux vector splitting approach is also applied at irregular grid points, 
which can be employed to stabilize the irregular grid stencils. Assuming subsonic flow speeds in the vicinity 
of the immersed boundary (no slip wall), upwinding can always be applied in the upwind direction towards 
the wall. This is in contrast to the scalar advection equation where the advection speed determines the 
upwind direction. The basic ideas of the IIM implementation are discussed for the 1D Euler equations in 
conservative variable formulation, 


p pu 
pu | =| pur+p |, (5) 
Be fs u( E+ p) 


where the total energy can be expressed as E = p (T'/((y — 1)yM?) + u?/2). The fluid is assumed to be an 
ideal gas, which can be modeled by the equation of state. Applying the flux vector splitting approach to the 
1D Euler equations leads to 

W, =F, =F, + Fr, (6) 


where the solution vector can be written as W = (w1,w2,w3)’ = (p,pu,E)’. In order to numerically 
analyze the 1D Euler system in Equation (6), the Jacobian, OF /OW, is introduced, 


OF OW OF” OW OFT OW 


Wi= aw ae = OW Oe + OW Ox” _ 
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This leads to an equation including the derivatives of the forward/backward-flux vectors, b=: with respect 
to the solution vector W. The gradient of the convective flux vectors, F*, is discretized with the discrete 
forward difference operator, D*, (for the forward-flux) and in backward direction with the discrete backward 
difference operator, D.. , (for the backward-flux), 


OF* OW OF OW 


“i= ow On 1 OW On is 
~ OF™ (pt.w) + 2F_. (pew) 
OW —~* OW ~*~ 


difference operators, De, are provided in Brehm et al.?! For an irregular grid point with an immersed 
boundary located to the left, the optimization procedure is only applied to the differentiation of the forward 
flux and to the differentiation of the backward flux when the immersed boundary located to the right. ‘The 
n'”-order accurate finite difference operator for the differentiation of the backward flux is obtained by utilizing 
n+ 2 grid points (including the boundary point) and applying a least-squares procedure to determine the 
stencil coefficients for the over-determined system of equations. In order to obtain pressure and temperature 
boundary conditions a weighted-least-squares interpolation approach is utilized. More details about the 


basics of the immersed boundary method can be found in Brehm et al.?°>?! 


III. Simulation Setup 





Figure 1. Un-installed open rotor configuration used in the present CFD analysis. 


Figure 1 presents the geometry used in the present CFD simulation. The CROR blade set shown in figure 
1 is a relatively modern GE design, namely F31/A31, with 12 blades on the front rotor and 10 blades on 
the second rotor. The simulation data is compared to the high-speed aerodynamics and acoustic data that 
was obtained for a model scale geometry in NASA’s 8-foot x 6-foot wind tunnel at cruise conditions. The 
front and aft rotor diameters are 0.66m and 0.63m, respectively. The two rotors are separated by 0.2m. The 
aerodynamic and acoustic data were acquired for the rotors in un-installed configuration without fuselage 
and pylon at zero angle of attack. 

A cut through the computational mesh employed for the current simulations is shown in Figures 2. 
Three grid resolutions are employed with 150, 250, and 500 million grid points to test the sensitivity of the 
simulation results with respect to the grid resolution. Note that all results shown in this abstract are based 
on the coarse grid simulations. Each Cartesian block (marked with solid black lines) contains 16° grid points. 
Grid refinement is applied in regions where high unsteadiness and large gradients are expected. Since storing 
the entire volume data for each time-step requires excessively large memory, several sampling planes were 
used to gather the unsteady flow data at each time-step. Most of the unsteady flow analysis was applied to the 
data recorded in these sampling planes. Note that throughout this paper all flow quantities were normalized 
with the free stream flow quantities, i.e., free-stream velocity, V,., pressure, P.., and temperature, T,,. The 
Reynolds number based on blade diameter D, V., and T,, is approximately 1.2 x 10". The free-stream 
Mach number is 0.78. 
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Figure 2. Computational mesh layout in the vicinity of the open rotor. 
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Figure 3. Colored gauge-pressure in the center plane and iso-surface of @-criterion. The black dots mark the Kulite 
sensor locations where experimental data is available. 










IV. Main Flow Features 


Color contours of gauge pressure and an iso-surface of Q-criterion is shown in figure 3. The black dots 
mark the Kulite sensor locations where experimental data is available. The Kulite sensors are numbered from 
left-to-right and top-to-bottom in increasing order. On the front and aft rotors, tip-vortices are generated at 
the end of each blade. The tip-vortices generated at the front blade impinge on the aft rotor blade as well as 
interact with the tip-vortices of the second rotor. The color contours of gauge pressure nicely visualize strong 
pressure waves that are generated by each rotor blade. These distinct pressure waves will generate strong 
peaks in the frequency spectra at distinct blade-passing-frequencies (BPF). One of the key challenges for 
simulating this flow field is to also capture broadband noise that is generated by the finer-scale turbulence in 
the wake. For lower rotor speeds, the broadband noise component becomes more important in comparison 
to the tonal noise and, thus, an accurate prediction is required. In order to capture the broadband noise, 
sufficient grid resolution is required to resolve the wide range of temporal and spatial scales in the turbulent 
wake. 
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Figure 4. Sound pressure level (SPL) spectra at Kulites (a) #9 and (b) #10 located at a vertical distance of 0.43m. 
Experimental data is shown with red lines and coarse grid computations with black lines. 


V. Preliminary Comparison with Experimental Data 


Figures 4a and 4b show a preliminary comparison of the coarse grid CFD results with experimental data 
obtained at different Kulite locations. It is expected that the coarse-grid resolution is not sufficient to capture 
some of the higher-order interactions. Due to insufficient grid resolution the large tip vortices cannot break 
up as quickly and, therefore, some of the distinct tones are slightly over-predicted in the spectra. This result 
is consistent with the observations by Housman et al. utilizing a conventional overset approach to simulate 
the flow around the open rotor. 


VI. Acoustic Analysis 


The following discussion will provide a brief overview of the acoustic analysis tools that will be used in the 
final version of the paper to identify the dominant noise sources and study the noise generation mechanisms. 
Some preliminary results are provided. 


VI.A. Proper Orthogonal Decomposition 


POD decomposes the flow field into a set of orthogonal basis functions that capture most of the flow energy 
in terms of a user-defined norm with the least number of modes.”° In many fluid dynamics applications, 
POD is used to identify energetic flow structures. A dynamical significance of existing POD modes is not 
guaranteed. Freund and Colonius?? applied POD to a Mach 0.9 turbulent jet. and noted that the sound 
radiated by turbulence amounts to a very small fraction of the total energy. ‘Thus, a POD approach in the 
temporal domain may not be able to capture the physics of the radiated sound waves. When analyzing 
the acoustic field it may be advantageous to follow an approach similar to the one described in Suzuki et 
al.2° Instead of using snapshots in time, Suzuki et al.2° used Fourier amplitude distributions as snapshots 
for POD. In the current study, we are dealing with mainly broadband noise sources as well as intermittent 
events and, therefore, POD was applied in the temporal domain. The frequency information of dominant 
unsteady jet and shock dynamics is obtained by applying a Fourier transform to the POD time-coefficients. 

Furthermore, the appropriate norm for capturing particular flow physics is not necessarily clear. Freund 
and Colonius?? have found that the effectiveness of a given representation of the initial flow data depends 
strongly upon the norm used and the data represented. For problems where the spatial dimension, m, is 
much greater than the number of available time-steps, NV, it is more efficient to employ the “snapshot” 
method?! for POD. Hence, in the current analysis the snapshot method is used. For more details about 
the POD snapshot method used here, see Chatterlee** for a basic introduction to POD, and Rowley** and 
Freund and Colonius?? on details on how it can be applied to fluid dynamics problems. In the current 
discussion, it is assumed that a vector of flow quantities, here q(#,t) = (p’,u’,v’, w’,s’), is defined over a 
region of interest 0 and additional weighting factors W weigh the contribution of the components of q(,t) 


7 of 15 


American Institute of Aeronautics and Astronautics 


in the vector norm, 
Ng 
a= [| Soonat |) ae, (9) 
2 \ e=1 


where N, is the number of components of g. By employing POD the series expansion converges optimally 
with respect to the L2-norm defined in Equation (9). Kinetic disturbance energy, the entropy fluctuation 
term, pressure or other combination of g can be used in the La-norm by choosing the weighting vector as 
3 = (0,1,1,1,0), (0,0,0,0,1), and (1,0,0,0,0), respectively. In Brehm et al.,3* it was shown that computing 
the most energetic p’-based POD modes is useful for obtaining a general sense of the unsteady shock motion 
for jet impingement problems. In order to study the noise generation process s’-based and K’-based POD 
modes provide some idea of how noise is generated in the source region when considering that the general 
structure of Lighthill’s stress tensor is composed of puju; and the entropy fluctuation term s’. More details 
about the computation of the POD modes and singular values associated with each mode are provided in 
appendix A. 


VIL.A.1. Preluminary POD Results 


Figures 5 display the most energetic pressure-based POD mode shapes (a,c,e,f) and corresponding spectra of 
the time-coefficients (b,d,f,h). Note that the POD analysis was conducted in a coordinate frame rotating with 
the front rotor. The tip vortex can be viewed as a steady mode in the rotating frame and is most dominant 
in the first two POD modes. The peak at shaft order 20 is associated with pressure wave generated by the aft 
rotor interacting with the front rotor and its tip vortices. In a fixed coordinate frame the aft blade generates 
pressure waves at shaft orders n x 10, where n is a positive integer value. In the contra-rotating frame (with 
the same rotational speed), the shaft orders are shifted to higher frequencies by a factor of 2. Considering 
the mathematical foundation of the POD modes it is clear that a single POD mode is essentially spatially 
stationary. This means that at least two POD modes are required in order to mathematically describe a 
traveling wave in terms of POD modes. This is the reason why POD modes frequently occur in pairs. The 
third and fourth most dominant POD modes display the aforementioned property. The frequency spectra 
of the time-coefficients are almost identical and the POD mode pair physically describes the propagation of 
high-amplitude pressure waves. 


VI.B. Spectral Near-Field Analysis 


Figure 7 shows the most dominant blade-passing-frequencies (BPF’) in the center-plane through the open 
rotor. Note that the results are expected to be significantly affected by the grid resolution. Hence, only 
preliminary conclusion can be drawn from the coarse grid simulation results. As one may expect, the most 
dominant pressure waves right in front of the open rotor are associated with the front rotor BPF1 while the aft 
rotor BPF2 remains dominant over a large region emerging from the back of the open rotor. Interestingly, a 
higher-order interaction (BPF1+2xBPF2) is dominant in the wake. The unsteady tail shock located slightly 
downstream of the abrupt change of the hub geometry generates significant pressure waves and appears to 
be dominated by BPF1. 

The noise generation process is further analyzed by considering the amplitude and phase plots at distinct 
frequencies. ‘The amplitude distribution for BPF1 and BPF2 display large amplitudes at the front and aft 
blade set, respectively. The tip vortices from the front and aft blades remain dominant in the wake and an 
interaction with the tail shock can be noted. Another weaker shock is formed right behind the aft blade and 
it appears to be dominated by BPF2. The amplitude and phase plots for BPF1+2xBPF2 seem to highlight 
dominant pressure waves originating from a region right behind the aft rotor in front of the aforementioned 
weaker shock. 
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Figure 5. Most energetic pressure-based POD mode shapes (a,c,e,g) and corresponding spectra of the time-coefficients 
(b,d,f,h) in a rotating frame with the front rotor. 
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Figure 6. Most dominant frequencies (in Hz) and corresponding Fourier amplitude displayed in the center-plane cutting 
through the open rotor. Note that the blade-passing-frequency of the front rotor is BPF1=1370Hz and BPF2=1140Hz. 
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Figure 7. Amplitude and phase distributions at different blade-passing-frequencies: (a) BPF1, (b) BPF2, and (c) 
BPF1+2xBPF2. 
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VI.C. Causality Method 


One of the key objectives of this work is to identify the dominant noise sources and uncover the relevant 
physical mechanisms responsible for noise generation. ‘To investigate noise generation mechanisms, the 
dynamics of the unsteady flow field discussed in section IV need to be connected to the radiated sound 
field. Simultaneous visualizations of the instantaneous flow and sound fields can provide an idea on specific 
events that generate noise, see for example.?° °’ In the current study, we employ a more direct approach by 
using the causality method, which relies on computing cross-correlations between flow quantities inside the 
open rotor flow field and the radiated acoustic pressure outside the non-linear flow regime. An attractive 
feature of this method is that the effects of scattering, absorption and refraction on sound radiation are 
automatically included by simultaneously extracting information from both the flow and acoustic fields. The 
present work is based upon the causality idea proposed in earlier works by.?° 42 Computing the correlation 
between turbulent fluctuations inside the unsteady open rotor flow field and the radiated noise in the far field 
is the most direct and unambiguous way of identifying the relevant noise sources. ‘The current theoretical 
framework is based on the acoustic analogy employed in Lighthill’s equation (10). It is well-known that the 
decomposition proposed in Lighthill’s analogy is non-unique. Flow-acoustic interactions are generally not 
distinguished from true acoustic sources.*! However, more elaborate formulations aimed at extending the 
acoustic analogy to account for the flow-acoustic interactions sacrifice the inherent simplicity of Lighthill’s 
original approach. In the current paper, flow-acoustic interaction is not explicitly accounted for and the 
noise sources shall be interpreted in terms of Lighthill’s sources. ‘These source term, however, contains fluid 
dynamics that do not follow the homogeneous-medium second-order wave equation. Because flow-acoustic 
interactions have been demonstrated to be important in shear layers, the results need to be interpreted 
carefully.4?: 49 


Lighthill’s equation is defined as 
Op 202, _ 
ee 


oT,;, 


1 
OD,02; ( 0) 


where the right-hand-side term contains Lighthill’s stress tensor T;; = pu;u; + 6:;(p — a2,p) absent the 
effect of viscosity, u; is the unsteady velocity field, and p is the unsteady pressure field. In the current 
discussion, it is assumed that the viscous effects are negligible for the noise generation as well as feedback 
from the acoustic field to the source. It is important to remember that Lighthill’s equation is simply a 
reformulation of the mass and momentum equations and is valid for every solution that obeys these equations. 
Lighthill’s equation cannot distinguish between radiating and non-radiating disturbances. Theoretically, the 
relationship dg, < w/a can be used to separate the radiating and non-radiating parts, where w is the angular 
frequency and a is the spatial wave number. In the current approach, we do not utilize this relationship 
explicitly and pursue a slightly different route to identify the radiating part of the solution. We follow the 
idea that multiplying with the acoustic pressure fluctuations in the far field acts simply as a filter operation, 
because sampling points located in the far field only sense the radiated part of disturbances.*4 

The most relevant parts of the causality approach are presented below. The theoretical foundation of 
our noise source identification strategy is laid out in more detail in Brehm et al.?7!_ The right-hand side 
of equation (10) can be regarded as a forcing term to the wave equation. Acoustically, the source terms 
represent a distribution of quadrupoles embedded in an ambient medium whereby flow discontinuities are not 
explicitly contemplated. The desired solution at a far-field point is derived in terms of an integral equation 
after manipulating equation (10) and utilizing the time-domain free-space Green’s function. 
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Following the analysis in Proudman* as more recently revisited by Panda et a we introduce the scalar 
component of the stress tensor J). that represent longitudinal quadrupoles. ‘The scalar quantity is the 
fluctuating stress tensor component T;. in the direction of the source point to the far-field observer location, 
i.e., P= af —2,. The acoustic intensity at a far-field point is finally obtained by multiplying equation (11) 
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with the far-field acoustic pressure and computing the acoustic intensity, 
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where the correlation function Cr,» is evaluated after shifting the near field data T;,. by the propagation time 
or time delay 7. Equation (12) provides an elegant way of linking the near field dynamics to the far-field 
acoustic pressure fluctuations. Applying such a cross-correlation is extremely powerful since it naturally 
separates the hydrodynamic and acoustic fluctuations.4* To avoid the computation of the second derivatives 
in time, we apply a temporal Fourier transform to equation (12) and find 


PSD, Care f) es [ CSDr..p’ (Db. f dv, (13) 


2 
TASS 


where the auto-correlation function of p’ is transformed to the power spectral density PSD, and the cross- 
correlation function between T;. and p’ turns into the cross-spectral density CSDr_»’. 
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Figure 8. Normalized cross-correlation of acoustic pressure at different sample points marked with white filled circle 
(e) and pressure on the sampling plane. Lines of constant time-delay are shown in the lower half of these figures with 
white solid lines. 


This method has previously been applied to a jet impingement problem discussed in Brehm et al.?! 


and it was shown to be very useful in studying the noise generation process. Figure 8 displays preliminary 
results utilizing normalized cross-correlations of acoustic pressure at different sample points marked with 
white filled circles and pressure on the sampling plane. The cross-correlation function varies from 0 to 1. 
Note the analysis employing the causality method will be significantly extended utilizing different terms of 
Lighthill’s stress tensor, normalized and non-normalized correlations. An important aspect of understanding 
the noise generation process is to determine which of the source terms are most dominant. Furthermore, to 
analyze turbo-machinery noise we will also study the contributions from the surface integrals in equation 
(11) and compare this result to traditional lower-order noise prediction methods that are based on a similar 
acoustic analogy. This way we will be able to point to weaknesses of these methods and possibly suggest 
improvements. 
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VII. Summary and Outlook 


In the present paper, we employ a higher-order accurate immersed boundary method to simulate the 
highly complex flow field around an open rotor. The current abstract provides an overview of the different 
topics that will be covered in more detail in the final paper. Two key aspects will be addressed in this paper: 
(1) the challenges of simulating the open rotor with our higher-order immersed boundary method and (2) a 
detailed analysis of the most dominant noise generation mechanisms of the open rotor propulsion system. The 
first part of the paper deals with the numerical simulation approach and the peculiarities of the immersed 
boundary method. In order to gain confidence in our simulation approach we are providing a detailed 
comparison with experimental data. In the second part of the paper, several post-processing tools, such as 
Fourier transforms, POD, and the causality method, are used to obtain a detailed understanding of the noise 
generation process. Finally, we want to emphasize that the results presented in the current abstract are 
based on a coarse grid simulation. The final paper will include a grid resolution study considering additional 
medium and fine grid solutions. 


Appendix A: Computation of POD Modes 


The POD eigenmodes are computed as 


N 
PO) (B) = S~ raz, ti), (14) 
1=0 


where the subscript 7 indicates the time-step (i = 0 to N), q(#,t;) is the vector of flow quantities at time- 


step 7, and ni) are individual elements of 7”). The right eigenvectors, 7”), are solutions to the algebraic 


eigenvalue problem: 
CA = MAM) | (15) 


where the correlation tensor C is defined in index notation as 


1 


Ci = 7 | 


Nz, ti), Q(z, ty) é (16) 


Here, the brackets denote the inner product [¢,q| = Jo oe) weak ) dx (see Equation (9)). The eigenfunc- 


tions are normalized with their corresponding eigenvalues BO), per] = bnn/A\™, and the time coefficients 


are computed as 


1 


a(t) = xy |atz,0), 9a) . (17) 





The flow field can be reconstructed/recovered from the eigenmodes and time coefficients: 


GE, t) » > a (y™ (2). (18) 


The POD modes are orthogonal to each other and sorted by their respective energy norm (see Equation 
(9)) contents, whereby the drop—off in mode energy toward the higher mode numbers is typically significant. 
Therefore, a small number of modes, J, will suffice for capturing a large percentage of the flow’s energy 
defined in Equation (9). This is particularly the case when the flow dynamics are governed by large energetic 
structures of organized (or coherent) fluid motion (such as in pure two-dimensional simulations). 
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